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Abstract 

We investigate the use of piecewise linear systems, whose coefficient matrix is a 
piecewise constant function of the solution itself. Such systems arise, for example, 
from the numerical solution of linear complementarity problems and in the numer- 
ical solution of free-surface problems. In particular, we here study their application 
to the numerical solution of both the (linear) parabolic obstacle problem and the ob- 
stacle problem. We propose a class of effective semi-iterative Newton-type methods 
to find the exact solution of such piecewise linear systems. We prove that the semi- 
iterative Newton- type methods have a global monotonic convergence property, i.e., 
the iterates converge monotonically to the exact solution in a finite number of steps. 
Numerical examples are presented to demonstrate the effectiveness of the proposed 
methods. 
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1 Introduction 



Due to their importance in both theory and applications, since the sixties a 
lot of interest has been devoted in the literature to both the obstacle and the 
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parabolic obstacle problems (see, e.g., [4,13]), which are nonlinear differen- 
tial problems of elliptic or parabolic type, respectively. Since the beginning, 
their theoretical study has been developed within the more general context 
of variational inequalities; the interested reader can refer, e.g., to [15], where 
their mathematical-physical introduction is given together with many related 
abstract theoretical results. 

In this paper we are in particular concerned with the numerical solution of the 
linear obstacle and linear parabolic obstacle problems of second order, that 
is we assume that the differential operators involved in the inequalities are 
of second order and that they are linear with respect to the unknown func- 
tion (observe that, despite of their name, these differential problems are still 
nonlinear because they involve differential inequalities instead of equalities). 
Their equivalent formulations as complementarity problems [15] is explicitly 
used in this paper for their introduction. 

The first numerical schemes for the obstacle problem were based on finite 
element discretizations of such problems combined with projected relaxation 
methods used for the solution of the discrete problem [5]. However, the con- 
vergence rate of such iterative schemes depends on the mesh refinement, and 
the position of the free boundary of the coincidence set (i.e., the set where the 
solution of the problem coincides with the obstacle) is not taken into account. 
In order to make the convergence rate independent of the mesh refinement, 
several multigrid algorithms have been proposed in the literature (see, e.g., 
[6,18]). Another way proposed in the literature to solve the obstacle problem 
is based on the use of an active set strategy which can also be combined with 
the multigrid approach (see, e.g., [7,11]). The scheme for the solution of the 
discrete problem is iterative (observe however that for the linear case con- 
vergence is obtained in finitely many steps and it is monotone) and at each 
step it makes a problem linearization by specifying the active and the inactive 
part of the unknowns. The active set strategy is used to define an outer itera- 
tion and each inner iteration requires the solution of a reduced linear system 
which can be efficiently implemented by the multigrid approach (see, e.g., [7]) 
or by a preconditioned conjugate gradient method when the linear case with 
Dirichelet boundary conditions is considered [8]. Inexact semismooth Newton 
methods have been developed in [10]. A further alternative is proposed in 
[17], where a different iterative approach is considered for the linear obstacle 
problem. In this case an iterative approximation of the contact region is used 
(moving obstacle). 

Even if in this paper we do not deal explicitly with mesh adaptation, we 
consider it an important aspect to be developed in the future in conjunction 
with our schemes. Relating to the coupled active set and multilevel approach, 
mesh adaptation is considered for the linear elliptic case in [8], where some a 
posteriori error estimates are reported. 
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Concerning the parabolic case, the use of a regularization technique combined 
with a Lagrange multiplier approach is proposed in [9], where some numerical 
results are presented for the Black-Scholes model for American options. In 
particular, such results are obtained by using a second order (in time and 
space) finite difference discretization of each regularized problem which leads 
to the solution of nonsmooth nonlinear equations which are numerically solved 
by a semismooth Newton method. Some interesting numerical results related 
to the linear parabolic obstacle problem are also reported in [1] where an Euler 
implicit time scheme is combined with a finite element spatial discretization. 
In such paper some a posteriori error estimates are used in order to control 
the mesh refinement, both in space and in time. Again, the discrete problem is 
solved by using a semismooth Newton method (see, e.g., [9]). The moving mesh 
method, based again on a posteriori error estimates, is the strategy suggested 
in [12] for both improving the accuracy and reducing the computational cost of 
the finite element numerical approximations of the parabolic obstacle problem 
which, otherwise, may have a very poor efficiency However, no numerical 
result is given in that reference. 

In this paper, we shall consider the numerical modeling and solution of lin- 
ear obstacle problems by means of piecewise linear systems (PLS, hereafter), 
which have been recently introduced and investigated in [2,3], with application 
to flows in porous media. PLS are linear systems, whose coefficient matrix is 
a piecewise constant function of the solution itself. They can be used for the 
efficient modeling of a number of real-life problems. In particular, we here con- 
sider their application for the numerical solution of the linear classical obstacle 
problem and its parabolic version. The procedure proposed for the numerical 
solution of the associated discrete obstacle problems is an iteration having a 
monotonic finite convergence behaviour. For completeness and clarity reasons, 
in the paper we specify that such application of PLS can be also formulated 
as a special case of the dual-active set strategy introduced in [11] and there 
studied under the assumption that the coefficient matrix characterizing the 
discrete problem is an M-matrix. However, in our opinion the PLS formula- 
tion of such iteration is an interesting alternative because of its compactness 
and because it allows us to analyze the convergence features of the method 
under less restrictive hypotheses which have to be assumed when obstacle 
problems with Neumann boundary conditions are dealt with. 

The paper is organized as follows. In Section 2 we investigate the classical 
obstacle problem. Then, in Section 3 we consider its modeling via PLS, whose 
numerical solution is investigated in Sections 4. In Section 5 the parabolic 
obstacle problem is considered, whose solution turns out to be a particular 
instance of what stated in Section 4. Section 6 contains some numerical ex- 
amples dealing with both Dirichlet and Neumann boundary conditions and, 
finally, Section 7 contains a few conclusions. 



3 



2 The (classical) obstacle problem 



We here consider the following special linear systems which involve nonsmooth 
functions of the solution itself, 

min{0,a?} + Tmax{0,a;} = b, (1) 



where x = (xj), b = (bi) G R n , with b a known vector, 



^ max{0, Xi} ^ 



max{0, a;} 



max{0, x n } 



^ min{0, Xi} ^ 



min{0, x} 



min{0,x n } 



and T G IR nxn is a (known) irreducible matrix, satisfying either one of the 
following properties: 

Tl: T is an M-matrix (i.e., it can be written as T = al — B, with B > O 
and p(B) < a), or 

T2: null(T T ) = span(v), null(T) = span(iu), with v,w > (compo- 
nentwise), and T + D is an M-matrix for all diagonal matrices D ^ O (i.e., 
D > O and D + O). 

Note that, if ^ = (&) G M. n is a given known vector, upon a suitable variable 
transformation, the following problems can be taken back to problem (1), 



mm{£,x} + Tirmx{£,x} = b, (2) 
max{£, x} + T min{£, a;} = b. (3) 



One important motivation, for solving problem (1), stands in the efficient 
numerical modeling of the the linear obstacle problem. In more details, let 
us consider the problem in its simplest form (see, e.g., [15] for more general 
formulations): 

-Au>f, u>i/j, (u- V>)(A« + /) = 0, in Q, (4) 



with suitable prescribed boundary conditions on dQ, where / is a known 
function and ip is the obstacle. 
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After a suitable finite difference discretization of problem (4), one obtains a 
corresponding discrete complementarity problem in the form 

Tu>f, u>V, (u-^) T (Tu- f) = 0, (5) 

where u is the unknown solution, / depends on the function / and on the 
boundary conditions, if) is the discrete representation of the obstacle, and T 
is a matrix satisfying either Tl, if u is specified in at least one point of dQ, 
or T2, otherwise, the previous problem can then be reformulated as 

Ty>b, y>0, y T (Ty-b) = 0, (6) 

where b = f — Tip. The following result then holds true. 

Theorem 1 If x is a solution of PLS (1), then y = max{0,a;} is a solution 
of (6). 

Proof Let s be a solution of (1). Clearly, max{0,a;} always satisfies the 
second inequality in (6). Then, concerning the first inequality and the com- 
plementarity condition, the following cases can occur, when considering the 
generic ith entry of x: 

• Xi < 0. Consequently, min{0, Xj} = Xj and max{0, Xj} = 0. Moreover, one 
has that the ith component of the first inequality in (6) is satisfied. In fact, 
by setting the zth unit vector: 

efTmaxjO, x} > min{0, + ej Tmax{0, x} = b^. 

• Xi > 0. In such a case, min{0,Xj} = and max{0,Xj} = Xj. Moreover, the 
ith component of the first inequality in (6) turns out to be an equality. In 
fact: 

efTmax{0, x} = min{0, x,} + ej Tmax{0, x} = 6j. 

Consequently, one concludes that y = max{0, x} satisfies all the inequalities 
in (6), as well as the complementarity condition. □ 



3 Modeling through PLS 

We here show how the nonsmooth equation (1) can be efficiently reformulated 
by means of a suitable PLS. In more details, for a given vector x e M. n , let 
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define the following diagonal matrix: 



P(x) 



p(xi) 



\ 



p(x ri 



\ 



J 



with p(xi 



1 if n > 0, 



otherwise. 



(7) 



The following straightforward results then hold true. 

Lemma 2 P(x)x = max{0, x}, [I — P(x)] x — min{0, x}. 

Lemma 3 System (1) is equivalent to the following PLS: 

[I - P(x)+TP(x)]x = b. (8) 

For sake of completeness, we also mention that problems (2)-(3) can be re- 
spectively reformulated as 



[/ - P ( {x) + TPs(x)} (x — £) = b—(I + T)£, 
[Pt(x) + T(I - Ps(x))} ( x -t) = b-(I+ T)t, 



(9) 
(10) 



where 



Pf(x) 



p(x 1 ) 



with p(xi) 



1 if Xi > & 



otherwise. 



4 The Newton-type iteration for the obstacle problem 



Some preliminary results are stated at first in order to derive a Newton-type 
procedure for solving the PLS (8) and prove its convergence. Their proof is 
straightforward and is, therefore, omitted. 

Lemma 4 Let T satisfy Tl. Then, for any diagonal matrix P, O < P < I, 
both matrices I — P + TP and I — P + PT are M -matrices and, therefore, 
(I - P + TPy 1 >0, (I - P + PTY 1 > O. Moreover, if in addition P ^ I, 
the same result continues to hold when T satisfies T2. 

It is to be noted that the left-hand side of system (8) is not everywhere dif- 
ferentiable. Nevertheless, a Newton-type method for solving system (8) can be 
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deduced, 



x k+i = x k _ fj _pk + Tpk 



fc = 0,1 



? J * * * 5 



I — P k + TP k ) x k - b\ , 
where the upper index denotes the iteration step and (see (7)) 

P° = 0, P k = P{x k ), k = 1,2, (12) 

This simplifies to the following Picard iteration, 

P° = 0, (i - P k + TP k )x k+1 = b, k = 0,1, (13) 



The following result provides a straightforward stopping criterion for the iter- 
ation. 

Lemma 5 If, for some k > 0, one gets 

(pfc+i _ p k ) x k + l = o, (14) 

then x* = x k+1 is an exact solution of problem (8). 
Proof Since (P fc+1 - P k )x k+1 = 0, one has 

(/ - P k + TP k ) x k+1 = (i - P k+1 + TP k+1 ) x k+1 = b, 
i.e., x k+1 solves (8). □ 

Remark 1 Actually, iteration (13) combined with the stopping criterion (14) 
can be formulated as a special case of the dual-active set strategy described 
by Algorithm Al in [11]. However, the next theorem shows that the compact 
matrix formulation here considered allows us a corresponding compact analysis 
of its convergence behaviour which is here extended to the case where T satisfies 
property T2 instead o/Tl. In addition, we observe that the PLS formulation 
(8) of the linear discrete obstacle problem allows us to get significant results 
about the existence and uniqueness of the solution of the discrete problem even 
for the extended case (see Theorem 12). 

The iteration is well-defined under the following conditions. 

Theorem 6 Let matrix T in system (23) satisfy Tl. Then, the matrix 

J pk _|_ rppk^ 
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is an M-matrix, and the iteration (13) is well defined for all k until conver- 
gence. IfT satisfies T2, the same result holds true, provided that 

v T b < 0. (15) 



Proof The thesis easily follows from Lemma 4, if we are able to prove that, 
when T satisfies T2 and P k ^ I, then: 

• either the exit condition (14) holds true, 

• or P k+1 zfz / ) so that the next iteration is well-defined. 

In the first case, by virtue of Lemma 5, x* = x k+1 is a solution of the problem, 
so that no further iterations are needed. In the second case, we observe that, 
from the definition (7), one readily shows that 

(pfc+i _ p^x^ 1 > o. 

If P k+1 = J ; then this would imply that, from (15) and (13), and considering 
that v > 0, 

< v T (P k+1 -P k )x k+1 = v T {I-P k )x k+l = v T {I-P k +TP k )x k+l = v T b < 0. 

Consequently, the exit condition (14) holds true, so that x k+1 is solution of 
problem (1). □ 

Next, we prove that the iteration (13) satisfies an important property of 
monotony. Before that, we state the following preliminary result, whose proof 
is straightforward and is, therefore, omitted. 

Lemma 7 By setting as usual x k = (x k ) and x k+1 = (x k+1 ), for k > 1 one 
has: 

(P k x k+1 > P k - 1 x k > 0) (x k > x k+1 > 0) =>> (P fc+1 >P k >0). 
Theorem 8 Let the hypotheses of Theorem 6 hold true. Then, 

P k+1 >P k >0, k = 0,1, (16) 



Proof For k = (16) trivially holds true, since P° = O. For k > 1, let us 
prove, according to Lemma 7, that 

P k x k + 1 > pk~l x k > Q 



Since, from (13), 

(I — P k + TP k )x k+l = {I- P^ 1 + TP k ^)x k = b, 
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one then obtains: 



(/ - P k + P k T) P k x k+1 

= P k (l - P k + TP k ) x k+1 = P k (i - P k ' x + TP k ~ l ) x k 
= (l-P k + P k T) P k ^x k + (P k - P k ~ 1 )x k . 

By considering that (I-P k + P k Ty 1 > O, ( y P k -P k ~ 1 )x k > 0, and P°x x = 0, 
(17) then follows. □ 

This result, allows to state the finite convergence of iteration (13). 

Corollary 9 Let T satisfy either Tl or T2. If T satisfies T2, assume that 
also (15) is satisfied. Then, iteration (13) converges in at most n steps. 

Proof The finite convergence easily follows from the fact that / > P k+1 > 
P k > O, and from the fact that, as soon as P k+1 = then the exit condition 
(14) is satisfied. Obviously, by considering that P° = O, this will happen in 
at most n steps. □ 

Remark 10 Even though Corollary 9 establishes the finite convergence of it- 
eration (13), nevertheless the corresponding upper bound may be large, when 
the dimension of the system is large. However, several numerical tests have 
shown that convergence can occur in just a few iterates (see, e.g., the numerical 
tests in Section 6). 

Next, we present a conclusion on the existence of a solution for problem (8). 
We need the following preliminary result. 

Lemma 11 With reference to matrix P defined in (7), for any two vectors 
x = (xi) and y = (yi), there exists a diagonal matrix, 

W = diag(wi,...,w„), 0<W<I, 

depending on x and y, such that 

P(x)x-P(y)y = W-(x-y). (18) 
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Proof For each i = 1,2, ... ,n, it follows from (7) that either one of the 
following four cases occurs: 

Xi, yi > p(xi) = p(yi) = 1 =^ Ui = \\ 

Xi,yi<0 => = = w< = 0; 

(19) 

^ > > y l p(xi) = 1, p(yi) =0 =► < Wj = < 1; 

Xi < < p(xi) = 0, p(yi) = 1 =>• < u { = < I. 

This proves the validity of (18). □ 
We can now state the following result. 

Theorem 12 Let T satisfy Tl. Then a solution to problem (8) exists and is 
unique. In the case where T satisfies T2, then a solution of problem (8): 

• exists and is unique when v T b < 0; 

• exists but is not unique when v T b = 0; 

• doesn't exist when v T b > 0. 

Proof Concerning the existence of a solution, the thesis follows from Corol- 
lary 9, when T satisfies Tl, or T satisfies T2 and (15) holds true. It remains 
to prove that no solution exists when T satisfies T2 and v T b > 0. Indeed, 
if such a vector x would exist, by considering that v > and taking into 
account (8), then 

< v T b = v T [I - P(x) + TP(x)]x = v T [I - P(x)}x < 0, 

which is clearly impossible. 

Concerning uniqueness, let x and y be two solutions of (8). Then 

b = [I - P{x) + TP(x)]x =[I- P{y) + TP{y)]y. 

By virtue of Lemma 11, this implies that 

M • (x - y) = (I - W + TW){x -y) = 0, 

for a suitable diagonal matrix W, O < W < I. Consequently, by taking into 
account Lemma 4: 

• if T satisfies Tl, then M is nonsingular and uniqueness (x = y) follows; 

• if T satisfies T2, then M is nonsingular, and uniqueness (x = y) follows, if 
and only if W ^ I. By considering the possible cases (19), this is equivalent 
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to requiring that at least one entry of one the two vectors x and y is negative. 
This is indeed the case, when v T b < 0, since 

> v T b = v T [I - P(x) + TP(x))x = v T [I - P(x)]x, 

which, by considering that v > and (see Lemma 2) [I — P{x)\x < 0, 
implies that at least one entry of x is negative. On the other hand, when 
v T b = 0, then 

= v T b = v T [I - P(x)]x = v T [I - P(y)]y, 

which implies that P(x) = P(y) = I. Consequently, one obtains Tx = Ty, 
i.e., 

x — y G null(T) = span(tu). 
In more details, if x = (xj) is a solution of (8) such that 

minxj = 0, 

i 

(observe that, since w > 0, such a solution always exists), then all solutions 
of (8) are given by 

x(a) = x + aw, a > 0. □ 

Remark 13 It is remarkable to observe that iteration (13) converges to a 
solution of problem (8) under the same hypotheses that guarantee its exis- 
tence. Moreover, convergence to a solution is guaranteed also when there is no 
uniqueness (i.e., when matrix T satisfies T2 and v T b — 0). 

For completeness, we mention that for solving problem (2), i.e. for solving the 
PLS (9), the corresponding iteration is: 

/f = 0, (/- P£ + TP?) (x k+1 - £) = b-(I + T)£, k = 0,1,..., 

where, see (11), Pj? = P^(x k ). The same iteration can be used for solving 
problem (3), i.e. for solving the PLS (10), if P k is replaced with P k = (I — P k ). 



5 The parabolic obstacle problem 

We now consider the application of PLS for the numerical solution of special 
linear systems, involving nonsmooth functions of the solution itself, in the 
form 

x + Tm&x{0,x} = b, (20) 
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where, as before, matrix T satisfies either Tl or T2. One important moti- 
vation, for solving problem (20), stands in the efficient numerical modeling 
of the linear parabolic obstacle problem. In more details, let us consider the 
problem in its simplest form (see, e.g., [14] for more general formulations): 

u t > Au + f, u > if), 

(21) 

(u — if))(ut — Au — f)=0, in fi, for t > 0, 

with suitable prescribed initial and boundary conditions at t — and on 
dQ. Here u t is the partial time derivative of the (unknown) solution u, f is 
a known function, and ijj is the (known) function describing the obstacle. A 
suitable implicit, finite difference discretization of problem (21), generates a 
corresponding discrete complementarity problem in the form 

u n+l + Tu n+l >u n + f, u n+1 > tf>, 

(u n+1 - ip) T {u n+1 + Tu n+l - u n - f) = 0, 

where u n is the discrete approximation at the nth time step (u° being specified 
by the initial condition), the vector / depends on the function /, on the 
boundary conditions and on the timestep, ij) is the discrete representation 
of the obstacle, and T is a matrix satisfying either Tl if, in the boundary 
conditions at the (n + l)st time step, u is specified in at least one point of 
dCl, or T2, otherwise PI By setting y = u n+1 — ij) and by defining a suitable 
(known) vector b, the previous problem, to be solved at each time step, can 
be reformulated as 

y + Ty>b, y>0, y T {y + Ty-b) = 0. (22) 
The following result then holds true. 

Theorem 14 If x is a solution of (20), then y = max{0,a?} is a solution of 
(22). 

Proof Let a; be a solution of (20). Clearly, max{0,£c} always satisfies the 
second inequality in (22). Then, concerning the first inequality and the com- 
plementarity condition, the following cases can occur, when considering the 
generic zth entry of x: 



As in the previous case, the problems in the literature generally prescribe the 
value of the solution at the boundary. For completeness, however, also in this case 
we consider the more general kind of boundary conditions. 
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• Xi < 0. Consequently, max{0,£j} = 0. Moreover, one has that the ith 
component of the first inequality in (22) is satisfied. Indeed, by setting 
the ith unit vector: 

max{0,:rj} + ejT max{0, x} > Xi + efTmax{0, x} = hi. 

• Xi > 0. In such a case, max{0, x{\ = x^. Moreover, the ith component of the 
first inequality in (22) turns out to be an equality. In fact: 

max{0,:Tj} + e^TmaxjO, x} = x% + efTmax{0, a?} = hi. 

One then concludes that y = max{0,a?} satisfies all the inequalities in (22), 
as well as the complementarity condition. □ 

From Lemma 2, the following straightforward result follows. 

Lemma 15 System (20) is equivalent to the following PLS: 

[I + TP(x)]x = b. (23) 



We now present an iterative procedure for solving (23), whose analysis will be 
taken back to what stated in Section 4. By using similar arguments to those 
used in that section, the iteration for solving (23) is given by: 

F° = 0, (i + TP k ) x k+1 = b, A; = 0,1,..., (24) 



where P k is given, as usual, by (12) and matrix T satisfies either Tl of T2. 
However, we observe that iteration (24) can be formally rewritten as 



P° = 0, I - P k + (I + T)P l 



x 



k+l 



b. 



k = 0,1, 



This implies that, since matrix (/ + T) is obviously Tl, for the iteration (24) 
hold all the results stated in the previous Section 4 for the iteration (13), 
in the Tl case. Namely, the iteration (24) always converges to the (unique) 
solution of problem (23) in a finite number of steps. 



6 Numerical tests 



In this section, we report a few numerical results for the above iterative meth- 
ods. In particular, in Section 6.1 we test the iteration (13) for solving the 
classical obstacle problem, whereas in Section 6.2 we test the iteration (24) 
for the parabolic obstacle problem. As it seems more frequent in the litera- 
ture, we mainly consider problems with Dirichlet boundary conditions and, 



13 



as a consequence, discrete problems with matrices T satisfying property Tl. 
However, in order to present also some results related to the extended case 
of a matrix T satisfying T2, a variant with Neumann boundary conditions is 
also examined, for both problems presented in the elliptic case. 

We specify that in all the reported tests, the linear systems associated with 
our iteration, which are sparse and nonsymmetric, are solved by using QMR 
(which is a standard iterative solver for unsymmetric systems; see, e.g., [16]) 

6.1 The obstacle problem 

Let us consider the following problem, whose obstacle is a non-smooth "tent- 
shaped" function: 

—Am > 0, u(x, y) > min(l — 2 — \y\) = ip(x, y), 

(25) 

Au(u - V) = 0, (x, y)en=(-l,l)x (-2, 2), u \ m = ^ 

Problem (25) is discretized, by using the standard 5-points second-order dif- 
ference scheme, on a cartesian grid with stepsizes 

2Ax = Ay = ~~~ -~ , (26) 
y N+ 1 y J 

The resulting PLS has then dimension n = N 2 . Because of the Dirichlet 
boundary conditions, the matrix T of such PLS, see Section 2, turns out 
to satisfy Tl. Figure 1 shows the plot of the computed numerical solution, 
whereas in Table 1 we list the number of iterations of (13), K, required to get 
convergence, for various values of N. 

A variant of the above problem is also considered. In this case we assume homo- 
geneous Neumann boundary conditions and a non vanishing forcing function 
/ constantly equal to — lj 3 l The Neumann boundary conditions have been 
discretized by using the standard 3-points second order forward an backward 
difference scheme. After eliminating the boundary unknowns, the associated 
PLS has always dimension n = N 2 but in this case the corresponding matrix 

2 We mention this only for sake of completeness: actually, it is not intended to 
discuss here the efficient solution of such linear systems which, by the way, can be 
further improved by using a suitable preconditioning technique. 

3 The last change has been introduced in order to deal with a discrete problem 
admitting a unique solution. 
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y 



-2 -1 



x 



Fig. 1. Solution of problem (25). 
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Fig. 3. Solution of problem (27) with Dirichlet boundary conditions (28), C = —20. 

T satisfies T2 and v = w with all unit components. The forcing term / is such 
that the right-hand side vector b in (8) has negative sum and, thus, Theorem 
12 allows us to state that the discrete problem admits a unique solution. Fig- 
ure 2 shows the plot of the computed numerical solution, whereas in Table 1 
we list the number of iterations of (13), Ky, required to get convergence, for 
various values of N. 

The second test problem is the elastic-plastic torsion problem in [17], 
—Au > C, u(x, y) > — min(x, 1 — x, y, 1 — y) = i/)(x, y), 



where C < is a given constant, and with homogeneous 
Dirichlet boundary conditions, 



It is known that the larger |C|, the more difficult the problem. Also in this 
case, a standard discretization, on a cartesian grid with stepsizes 



(x,y)en=(0 1 l) 2 , 



(27) 



(Au + C){u-i/>) = 0, 



u \an= i> |sn= 0. 



(28) 



Aa; = Ay 



1 



(29) 



N+l 



1(3 



Fig. 4. Coincidence set (dotted region) related to problem (27) with Dirichlet bound- 
ary conditions (28). On the left for C = —5 and on the right for C = —20. 



Fig. 5. Coincidence set (dotted region) related to problem (27) with Neumann 
boundary conditions (30). On the left for C = —5 and on the right for C = —20. 



leads to a PLS of dimension n = N 2 , whose matrix T satisfies Tl. Figure 3 
shows the plot of the computed numerical solution (C = —20), whereas in 
Table 2 we list the number of iterations, K(C), required to get convergence 
for different values of C and of the discretization parameter iV in (29). In 
such a case, the number of the required iterations of (24) decreases, as \C\ 
increases. This behaviour can be explained considering that, as \C\ increases, 
the coincidence set, i.e. the set of points in Q where u — ip, enlarges, as shown 
in Figure 4. Consequently, our initialization becomes nearer to the solution. 
In fact, in our iterative procedure (13) we assume P° = O; thus x 1 = b and, 
for sufficiently fine grids, the negativity of C and the analytical expression of 
the obstacle function imply that all, or almost all, the components of b are 
negative. Thus, as the solution u of (5) is initialized with maxjO^ 1 } + if), 
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Table 1 

Numerical results for problem (25) and its Neumann variant, both discretized with 
stepsizes (26); K and Ky respectively denote the number of iterations of (13) for 
problem (25) and for its variant. 



N 


25 


50 


75 


100 


n 


625 


2500 


5625 


10000 


K 


6 


10 


10 


12 


Ky 


12 


25 


37 


49 



Table 2 

Numerical results for problem (27) discretized with stepsizes (29), either with 
Dirichlet boundary conditions (28) or with Neumann boundary conditions (30); 
K(C) denotes the number of iterations of (13) for the specified value of the param- 
eter C. 



N 


25 


50 


75 


100 


n 


625 


2500 


5625 


10000 


K(C = 


-5) 


9 


17 


25 


32 


K{C = 


-10) 


5 


10 


13 


16 


K{C = 


-15) 


4 


7 


9 


11 


K{C = 


-20) 


4 


5 


7 


9 



its initialization in this case is almost everywhere coincident with the obstacle 
and, thus, closer and closer to the final solution as \C\ is increased. 

Even for this problem we have considered a variant with Neumann boundary 
conditions, which have been chosen in order to get a solution with a shape 
analogous to that of the solution of the Dirichlet problem. In more detail, (28) 
is replaced by the following non homogeneous Neumann boundary conditions, 



du 
dn 



tm 



dip 
dn 



(30) 



where a suitable extension of the derivative of i/j at the corner points of the 
domain is used. By using again a second order discretization of the boundary 
conditions and eliminating the boundary unknowns, the problem dimension is 
unchanged and, as for the variant of the first problem, we deal with a matrix 
T satisfying T2, and v = w with all unit components. Also in this case, from 
Theorem 12, we obtain a unique solution for the associated discrete problem. 
For all the values of N and C considered in Table 2, as outlined in the caption 
of the table, convergence has been obtained with the same number of iterations 
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as for the Dirichlet case. No plot of the solution is reported for this problem 
with Neumann conditions because, for all the considered values of C, it is 
analogous to that related to the Dirichlet case. We only present the obtained 
coincidence sets in Figure 5 which, when compared with the corresponding 
ones in Figure 4, better show the small difference near the boundary between 
the obtained solution and that related to the Dirichlet case. 

6.2 The parabolic obstacle problem 

The problems that we shall consider here, are evolutionary versions of those 
considered in Section 6.1. In more details, the first problem is given by 

u t > Au, u(x, y, t) > ijj(x, y), 

(x,y,t) e Q x (0,r], u \ dn = -, 

where ip and Q are the same items defined in (25). The spatial discretization is 
the same used for that problem (see (26)), whereas the discretization in time 
is, for sake of simplicity, done by means of the implicit Euler method, by using 
a constant stepsize 

At = -, (32) 



(u t - Au){u-ip) = 0, 



u |t=o= max ip, 




being v the number of time steps. The resulting PLS, to be solved at each time 
step, has dimension n = N 2 , whose matrix T is the same as that obtained for 
problem (25). Table 3 summarizes the obtained results, in terms of required 
number of iterations, for r = 10 4 and v = 20. In such a case, the approximation 
at t = t is quite close to the limit solution plotted in Figure 1. Here the number 
of iterations of (24) appears to be independednt of spatial resolution. 

The second test problem is the evolutionary version of the elastic-plastic tor- 
sion problem (27): 

u t > Au + C, u(x, y,t) > ip(x,y), {u t — Au — C){u — iff) =0, 

(33) 

(x,y,t) e Q x (0,r], u \an= 0, u | t =o= max (ip, 0) = 0, 

where tj) and Q are the same items defined in (27). The spatial discretization is 
the same used for that problem (see (29)), whereas the discretization in time 
is, for sake of simplicity, done by means of the implicit Euler method, by using 
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Table 3 

Number of iterations for problem (31), discretized with stepsizes (26) and (32), at 
each timestep iAt, i = 1, . . . , 20. 



N 


25 


50 


75 


100 


i\n 


625 


2500 


5625 


10000 


1 


5 


6 


8 


8 


2 


5 


6 


6 


6 


20 


5 


6 


6 


6 



Table 4 



Number of iterations for problem (33), discretized with stepsizes (29) and (32), at 
each timestep iAt, i = 1, . . . , 20. 



N 


25 


50 


75 


100 


n 


625 


2500 


5625 


10000 


i\C 


-5 -10 -15 -20 


-5 -10 -15 -20 


-5 -10 -15 -20 


-5 -10 -15 -20 


1 
20 


9 5 4 4 
9 5 4 4 


17 10 7 5 
17 10 7 5 


25 13 9 7 
25 13 9 7 


32 16 11 9 
32 16 11 9 



a constant stepsize (32). Also in this case, the resulting PLS, to be solved at 
each time step, has dimension n = N 2 , whose matrix T is the same as that 
obtained for problem (27). Table 4 summarizes the obtained results, in terms 
of required number of iterations, for r = 5 and v = 20. In such a case, the 
approximation at t = r is quite close to the limit solution plotted in Figure 3. 
The number of iterations required for obtaining the solution turns out to be 
quite similar to that listed in Table 2 for the corresponding stationary problem. 



7 Conclusions 

Two simple semi- iterative Newton-type procedures for solving certain classes 
of piecewise linear systems have been investigated. Such piecewise linear sys- 
tems are derived from the efficient modeling of obstacle problems. It has been 
shown that, under rather general assumptions, the iterates are well defined 
and monotonically converge to an exact solution of the given system in a fi- 
nite number of steps. A few numerical examples, concerning both the classical 
obstacle problem, and its evolutionary (i.e., parabolic) counterpart, prove the 
effectiveness of the proposed methods. 
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